NATRE microstructure dataset#

Study NATRE dataset, reproduce Ferrari & Polzin (2005)

%load_ext watermark
%matplotlib inline

import glob

import cf_xarray as cfxr
import dcpy
import distributed
import gsw_xarray
import hvplot.xarray
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import tqdm
import xfilter
from IPython.display import Image

import eddydiff as ed
import xarray as xr

xr.set_options(keep_attrs=True)

plt.rcParams["figure.dpi"] = 140
plt.rcParams["savefig.dpi"] = 200
plt.style.use("ggplot")

%watermark -iv
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/statsmodels/compat/pandas.py:65: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
  from pandas import Int64Index as NumericIndex
numpy      : 1.22.3
tqdm       : 4.64.0
distributed: 2022.4.0
hvplot     : 0.7.3
cf_xarray  : 0.7.1.dev16+g7305fd8
gsw_xarray : 0.3.0
xarray     : 2022.6.0
matplotlib : 3.5.1
eddydiff   : 0.1
dcpy       : 0.2.1.dev11+gebbbf87.d20220428
xfilter    : 0.2.1.dev6+g5fdc006
if "client" in locals():
    client.cluster.close()
    client.close()
client = distributed.Client(n_workers=2, threads_per_worker=2)
client

Client

Client-8948807e-35dd-11ed-95d3-3af9d394f1c6

Connection method: Cluster object Cluster type: distributed.LocalCluster
Dashboard: http://127.0.0.1:8787/status

Cluster Info

Image("../images/natre-large-scale.png")
_images/natre-examine_3_0.png

Read dataset#

created using eddydiff.natre.combine_natre_files()

natre = ed.natre.read_natre(load=True)
13 40
0  values == 0
28090
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/numba/np/ufunc/gufunc.py:151: RuntimeWarning: invalid value encountered in get_gradient_sign
  return self.ufunc(*args, **kwargs)
ed.plot.histogram_turb_estimates(natre)

\(R_ρ, K_T, K_ρ\)#

(
    natre..sel(pres=slice(250, None))
    .where(np.abs(natre.) < 100)
    .mean(["latitude", "longitude"])
    .reset_coords(drop=True)
    .hvplot()
).opts(show_grid=True)
ds = xr.Dataset(
    {"KrhoKt": (natre.Krho / natre.Kt), "Rρ": natre..where(np.abs(natre.) < 5)}
).sel(pres=slice(250, None))
ds.plot.scatter(y="KrhoKt", x="Rρ")
plt.yscale("log")
_images/natre-examine_9_0.png
import colorcet
import datashader as dshdr
from holoviews.operation.datashader import datashade
df = ds.to_dataframe()
df["KrhoKt"] = np.log10(df.KrhoKt)
pts = hv.Points(df, ["KrhoKt", "Rρ"])
datashade(pts, cmap=colorcet.fire).opts(width=600, height=600)

TS plot#

hdl, axes = dcpy.oceans.TSplot(
    natre.salt[::2],
    natre.theta[::2],
    Pref=1000,
    # rho_levels=bins,
    Sbins=np.arange(35, 35.5, 0.025),
    Tbins=np.arange(4, 10, 0.3),
    hexbin=False,
    plot_kwargs={"alpha": 0.01},
)
_images/natre-examine_13_0.png
(
    natre["chi"].rolling(depth=100, center=True, min_periods=1).mean().compute()
).plot.line(
    y="depth",
    col="longitude",
    row="latitude",
    yincrease=False,
    xlim=[1e-12, 1e-7],
    ylim=(2000, 0),
    xscale="log",
)
tornado.application - ERROR - Exception in callback functools.partial(<bound method IOLoop._discard_future_result of <zmq.eventloop.ioloop.ZMQIOLoop object at 0x7f70c8187340>>, <Task finished name='Task-481' coro=<Cluster._sync_cluster_info() done, defined at /home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/deploy/cluster.py:104> exception=OSError('Timed out trying to connect to tcp://127.0.0.1:45051 after 3 s')>)
Traceback (most recent call last):
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/comm/core.py", line 284, in connect
    comm = await asyncio.wait_for(
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/asyncio/tasks.py", line 498, in wait_for
    raise exceptions.TimeoutError()
asyncio.exceptions.TimeoutError

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/tornado/ioloop.py", line 741, in _run_callback
    ret = callback()
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/tornado/ioloop.py", line 765, in _discard_future_result
    future.result()
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/deploy/cluster.py", line 105, in _sync_cluster_info
    await self.scheduler_comm.set_metadata(
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/core.py", line 785, in send_recv_from_rpc
    comm = await self.live_comm()
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/core.py", line 742, in live_comm
    comm = await connect(
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/comm/core.py", line 308, in connect
    raise OSError(
OSError: Timed out trying to connect to tcp://127.0.0.1:45051 after 3 s
tornado.application - ERROR - Exception in callback functools.partial(<bound method IOLoop._discard_future_result of <zmq.eventloop.ioloop.ZMQIOLoop object at 0x7f70c8187340>>, <Task finished name='Task-488' coro=<Cluster._sync_cluster_info() done, defined at /home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/deploy/cluster.py:104> exception=OSError('Timed out trying to connect to tcp://127.0.0.1:45051 after 3 s')>)
Traceback (most recent call last):
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/comm/core.py", line 284, in connect
    comm = await asyncio.wait_for(
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/asyncio/tasks.py", line 498, in wait_for
    raise exceptions.TimeoutError()
asyncio.exceptions.TimeoutError

The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/tornado/ioloop.py", line 741, in _run_callback
    ret = callback()
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/tornado/ioloop.py", line 765, in _discard_future_result
    future.result()
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/deploy/cluster.py", line 105, in _sync_cluster_info
    await self.scheduler_comm.set_metadata(
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/core.py", line 785, in send_recv_from_rpc
    comm = await self.live_comm()
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/core.py", line 742, in live_comm
    comm = await connect(
  File "/home/deepak/miniconda3/envs/dcpy/lib/python3.8/site-packages/distributed/comm/core.py", line 308, in connect
    raise OSError(
OSError: Timed out trying to connect to tcp://127.0.0.1:45051 after 3 s
<xarray.plot.facetgrid.FacetGrid at 0x7f70a01c4e50>
_images/natre-examine_14_2.png

Procedure#

The large-scale average operator ⟨⟩ represents

  1. a horizontal average over the survey lateral scale,

  2. a vertical average over O(100) m, and

  3. a time average over the 18-day survey.

The mean fields are derived by averaging all variables along neutral-density surfaces \(γ_n\)

Mean vertical gradients \(∂_z θ_m\), \(∂_z S_m\), and \(∂_z b_m\), are calculated from O(100)-m linear fits to \(θ_m(z_n)\), \(S_m(z_n)\), and \(b_m(z_n)\), where \(z_n\) is the mean depth of each surface \(γ_n\).

large-scale temperature gradient \(|∇_n θ_m|\) is estimated by fitting a plane to the 100 stations in the 400 km × 400 km large-scale survey grid and rms is computed as the standard deviation of the departures from the plane fit.

In theory, these bins are approx. 100m apart in neutral density

Choose bins#

# bins = (
#    natre.gamma_n.mean(["latitude", "longitude"])
#    .interp(depth=np.arange(150, 2001, 100))
#    .dropna("depth")
#    .data
# )

bins = ed.sections.choose_bins(natre.gamma_n, depth_range=np.arange(150, 2001, 100))
natre.gamma_n.stack({"latlon": ["latitude", "longitude"]}).drop("latlon").cf.plot(
    hue="latlon",
    y="depth",
    lw=0.5,
    add_legend=False,
    yincrease=False,
)
dcpy.plots.linex(bins)
_images/natre-examine_18_0.png

Make estimate#

Vertical gradients: ed.sections.fit1D

Mean vertical gradients \(∂_z θ_m\), \(∂_z S_m\), and \(∂_z b_m\), are calculated from O(100)-m linear fits to \(θ_m(z_n)\), \(S_m(z_n)\), and \(b_m(z_n)\), where \(z_n\) is the mean depth of each surface \(γ_n\).

I don’t understand this. The bins are O(100m) apart, so how do you “fit” straight lines over O(100)-m to a profile that has points every O(100)m.

Horizontal gradients: ed.sections.fit2D

large-scale temperature gradient \(|∇_n θ_m|\) is estimated by fitting a plane to the 100 stations in the 400 km × 400 km large-scale survey grid and rms is computed as the standard deviation of the departures from the plane fit.

Ke is totally wrong!, I need to fix my plane fitting

Bin-average in density bins#

micro = ed.sections.bin_average_vertical(
    natre.reset_coords("pres"), "neutral_density", bins
)
micro
<xarray.Dataset>
Dimensions:             (gamma_n: 17, bounds: 2)
Coordinates:
  * gamma_n             (gamma_n) float64 26.79 26.96 27.1 ... 27.93 27.95 27.96
    pres                (gamma_n) float64 303.1 407.2 ... 1.826e+03 1.902e+03
    reference_pressure  int64 1000
    num_obs             (gamma_n) int64 19763 19424 17093 ... 18916 21881 10032
    gamma_n_bounds      (bounds, gamma_n) float64 26.69 26.88 ... 27.96 27.97
Dimensions without coordinates: bounds
Data variables: (12/18)
    chi                 (gamma_n) float64 6.112e-09 3.64e-09 ... 2.707e-10
    eps                 (gamma_n) float64 6.149e-10 5.53e-10 ... 1.16e-10
    salt                (gamma_n) float64 36.16 35.87 35.68 ... 35.13 35.1 35.09
    temp                (gamma_n) float64 15.64 13.83 12.45 ... 4.776 4.47 4.271
    theta               (gamma_n) float64 15.75 13.92 12.52 ... 4.393 4.187
    pden                (gamma_n) float64 1.031e+03 1.031e+03 ... 1.032e+03
    ...                  ...
    Kt                  (gamma_n) float64 2.416e-05 2.335e-05 ... 1.511e-05
    KtTz                (gamma_n) float64 1.398e-07 1.077e-07 ... 3.357e-08
    dTdz_m              (gamma_n) float64 0.01861 0.01519 ... 0.002752 0.002534
    N2_m                (gamma_n) float64 1.95e-05 1.679e-05 ... 1.355e-06
    Krho_m              (gamma_n) float64 6.307e-06 6.589e-06 ... 1.713e-05
    Kt_m                (gamma_n) float64 8.828e-06 7.89e-06 ... 2.107e-05
Attributes: (12/13)
    Conventions:           CF-1.6
    netcdf_version:        4
    project:               North Atlantic Tracer Release Experiment (NATRE)
    expocode:              32OC250_4
    cast_number:           3.0
    title:                 Microstructure profiler data from the ship Oceanus...
    ...                    ...
    latitude:              27.533166666666666
    longitude:             -30.723333333333333
    chief_scientist:       Raymond W. Schmitt
    data_originator:       Polzin
    institution:           WHOI
    data_assembly_center:  CCHDO

Replicate some Ferrari & Polzin (2005) figures. \(K_ρ, ε\) look OK. \(K_e\) is wrong, something is wrong with my slope estimate / plane fitting

f, ax = plt.subplots(1, 3, sharey=True, constrained_layout=True)
micro.drop_vars("pres").Krho_m.cf.plot.step(xscale="log", xlim=(4e-6, 1e-4), ax=ax[0])
micro.drop_vars("pres").eps.cf.plot.step(ax=ax[1], xscale="log")
# chidens["Ke"].cf.plot(ax=ax[2])
[<matplotlib.lines.Line2D at 0x7f706d7c6e20>]
_images/natre-examine_26_1.png

Constructing mean profiles#

I experimented with constructing a mean profile using a large number of γ_n bins, and then fitting straight lines in approx 100m bins.

I have instead decided to do a similar fitting procedure in ed.sections.fit1D for each chosen O(100m) wide density bin. Otherwise I would need to figure out how to go from \(θ_z\) in this mean profile to an appropriate value for the density bin.

bins = ed.sections.choose_bins(
    natre.gamma_n.mean(["latitude", "longitude"]), np.arange(150, 2001, 10), decimals=6
)

props = natre[["temp", "salt", "gamma_n"]].interpolate_na("depth")
props["depth_"] = ("depth", props.depth.data)
mean_props = props.groupby_bins("gamma_n", bins).mean()
sliding = (
    mean_props  # .swap_dims({"gamma_n_bins": "depth_"})
    # .rename({"depth_": "depth"})
    .rolling(gamma_n_bins=10, center=True).construct("bin")
)
Δz = sliding.depth_.isel(bin=[-1, 0]).diff("bin").squeeze()
gradients = sliding.polyfit("bin", deg=1).sel(degree=1) / Δz
mean_props["Tz"] = gradients.temp_polyfit_coefficients
mean_props["Sz"] = gradients.salt_polyfit_coefficients
mean_props["N2"] = -9.81 / 1030 * gradients.gamma_n_polyfit_coefficients
mean_props

Mean gradients in a chosen O(100m) γ_n bin#

bins = ed.sections.choose_bins(natre.gamma_n, depth_range=np.arange(150, 2001, 100))
binned = natre.reset_coords().groupby_bins("gamma_n", bins=bins)

i = 0
for label, group in binned:
    i = i + 1
    if i == 5:
        break
ds = group.unstack()
ds
<xarray.Dataset>
Dimensions:             (latitude: 10, longitude: 10, depth: 547)
Coordinates:
  * latitude            (latitude) float64 27.5 27.1 26.7 ... 24.7 24.3 23.9
  * longitude           (longitude) float64 -30.7 -30.3 -29.8 ... -27.2 -26.8
  * depth               (depth) float64 594.6 595.0 595.8 ... 495.8 497.8 499.8
Data variables: (12/18)
    chi                 (latitude, longitude, depth) float64 1e-15 ... 6.867e-11
    eps                 (latitude, longitude, depth) float64 7.443e-11 ... 2....
    salt                (latitude, longitude, depth) float64 35.56 ... 35.62
    temp                (latitude, longitude, depth) float64 11.7 11.7 ... 11.89
    gamma_n             (latitude, longitude, depth) float64 27.16 ... 27.16
    pres                (latitude, longitude, depth) float64 599.5 ... 503.5
    ...                  ...
    chi_masked          (latitude, longitude, depth) float64 1e-15 ... 6.867e-11
    Krho                (latitude, longitude, depth) float64 1.744e-06 ... 4....
    KrhoTz              (latitude, longitude, depth) float64 1.276e-08 ... 6....
    eps_chi             (latitude, longitude, depth) float64 3.988e-16 ... 1....
    Kt                  (latitude, longitude, depth) float64 9.347e-12 ... 1....
    KtTz                (latitude, longitude, depth) float64 6.836e-14 ... 2....
Attributes: (12/13)
    Conventions:           CF-1.6
    netcdf_version:        4
    project:               North Atlantic Tracer Release Experiment (NATRE)
    expocode:              32OC250_4
    cast_number:           3.0
    title:                 Microstructure profiler data from the ship Oceanus...
    ...                    ...
    latitude:              27.533166666666666
    longitude:             -30.723333333333333
    chief_scientist:       Raymond W. Schmitt
    data_originator:       Polzin
    institution:           WHOI
    data_assembly_center:  CCHDO
ed.sections.average_density_bin(group).dTdz_m
-0.011270900780997465
0.012996880869341362
<xarray.DataArray 'dTdz_m' ()>
array(0.0112709)
Attributes:
    name:     $∂_z θ_m$
    units:    °C/m
group.unstack().depth
<xarray.DataArray 'depth' (depth: 547)>
array([594.6, 595. , 595.8, ..., 495.8, 497.8, 499.8])
Coordinates:
  * depth    (depth) float64 594.6 595.0 595.8 596.2 ... 493.4 495.8 497.8 499.8
ed.sections.fit1D(group, "theta", "depth", debug=True)
<xarray.DataArray 'depth' (depth: 10)>
array([551.590043, 561.407269, 572.258639, 583.53746 , 593.436723, 604.726001,
       616.072442, 627.213435, 637.901493, 648.193224])
Coordinates:
    gamma_n_bins  (depth) object (27.16, 27.173] ... (27.277, 27.29]
  * depth         (depth) float64 551.6 561.4 572.3 583.5 ... 627.2 637.9 648.2
-0.011354461927299291
array(-0.01135446)
_images/natre-examine_35_2.png

Computing along latitude and longitude lines#

def compute_single_line(section):
    return ed.sections.bin_average_vertical(
        section.cf.stack({"cast": ("latitude", "longitude")}).drop("cast"),
        "neutral_density",
        bins,
    )


lonlines = natre.groupby("longitude", squeeze=False).map(compute_single_line)
latlines = natre.groupby("latitude", squeeze=False).map(compute_single_line)

latlines.to_netcdf("../datasets/natre-along-lat-lines.nc")
lonlines.to_netcdf("../datasets/natre-along-lon-lines.nc")

Buoyancy Reynolds Number#

There are a bunch of low values

natre["ν"] = dcpy.oceans.visc(natre.salt, natre.temp, natre.pres)
natre["κ"] = dcpy.oceans.tdif(natre.salt, natre.temp, natre.pres)
natre["Reb"] = natre.eps / natre.ν / natre.N2
natre["Reb"].attrs["long_name"] = "$Re_b$"
del natre.Reb.attrs["units"]
natre["χmol"] = 2 * natre.κ * natre.Tz**2
_, bins, _ = np.log10(natre.χmol + natre.chi).plot.hist(
    bins=101, histtype="step", density=True, lw=2, aspect=3, size=3
)
np.log10(natre.chi).plot.hist(bins=bins, histtype="step", density=True, lw=2);
_images/natre-examine_41_0.png
kwargs = dict(bins=np.linspace(-2, 5, 101), density=True, histtype="step", lw=2)

np.log10(natre.Reb.sel(latitude=24.7, method="nearest")).plot.hist(
    **kwargs, aspect=3, size=3
)
np.log10(natre.Reb.sel(latitude=25.5, method="nearest")).plot.hist(**kwargs)

dcpy.plots.linex(np.log10([16, 200]))
/Users/dcherian/work/python/xarray/xarray/core/computation.py:734: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)
/Users/dcherian/work/python/xarray/xarray/core/computation.py:734: RuntimeWarning: invalid value encountered in log10
  result_data = func(*input_data)
_images/natre-examine_42_1.png

NATRE \(ε\) vs \(ε_χ\)#

np.log10(
    natre.eps_chi.where(np.abs(natre.Tz) > 5e-3).where(natre.eps_chi < 1e-6)
).plot.hist(bins=101, histtype="step", density=True)
np.log10(natre.eps).plot.hist(bins=101, histtype="step", density=True, yscale="linear");
_images/natre-examine_44_0.png
natre.eps_chi.where(np.abs(natre.Tz) > 5e-3).where(natre.eps_chi < 5e-6).mean(
    ["latitude", "longitude"]
).coarsen(pres=100, boundary="trim").mean().plot()
natre.eps.mean(["latitude", "longitude"]).coarsen(
    pres=100, boundary="trim"
).mean().plot(yscale="log")
[<matplotlib.lines.Line2D at 0x16f00e0a0>]
_images/natre-examine_45_1.png

Compute NATRE microscale budget terms#

In theory, these bins are approx. 100m apart in neutral density

natre = ed.natre.read_natre(load=True).drop("depth")
bins = ed.sections.choose_bins(natre.gamma_n, depth_range=np.arange(150, 2001, 100))
natre.gamma_n.stack({"latlon": ["latitude", "longitude"]}, create_index=False).cf.plot(
    hue="latlon",
    y="pres",
    lw=0.5,
    add_legend=False,
    yincrease=False,
)
dcpy.plots.linex(bins)
_images/natre-examine_47_0.png
micro.load(scheduler=client).to_netcdf("../datasets/natre-density-binned.nc")
OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
OMP: Info #271: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.

T profiles#

ds = natre.stack(
    {"station": ("latitude", "longitude")}, create_index=False
).assign_coords(station=np.arange(100))
import xfilter

Zname = ds.cf.axes["Z"][0]

CT = ds.CT.interpolate_na(Zname)
sortlo = xfilter.lowpass(CT, coord=Zname, freq=1 / 20)
sortlo100 = xfilter.lowpass(CT, coord=Zname, freq=1 / 100)

(
    sortlo100.hvplot.line(y=Zname, color="k")
    # * unsorted.CT.hvplot.line(y=Zname)
    * sortlo.hvplot.line(y=Zname, color="r")
    # * sort.CT.hvplot.line(y="pressure", ylim=(2000, 0))
).opts(ylim=(1200, 800), xlim=(6, 10), frame_width=500, frame_height=400)

Debugging estimates#

natre = ed.natre.read_natre(load=True, stack=False)
natre
13 40
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/numba/np/ufunc/gufunc.py:151: RuntimeWarning: invalid value encountered in get_gradient_sign
  return self.ufunc(*args, **kwargs)
0  values == 0
28090
<xarray.Dataset>
Dimensions:    (latitude: 10, longitude: 10, pres: 3981)
Coordinates:
  * latitude   (latitude) float64 27.5 27.1 26.7 26.3 ... 25.1 24.7 24.3 23.9
  * longitude  (longitude) float64 -30.7 -30.3 -29.8 -29.4 ... -27.6 -27.2 -26.8
  * pres       (pres) float64 10.0 10.5 11.0 11.5 ... 1.999e+03 2e+03 2e+03
    time       (latitude, longitude) datetime64[ns] 1992-03-28T15:28:59.99999...
    depth      (latitude, longitude, pres) float64 nan nan ... 1.978e+03
Data variables: (12/22)
    chi        (latitude, longitude, pres) float64 nan nan ... 1.145e-11
    eps        (latitude, longitude, pres) float64 nan nan ... 2.916e-11
    salt       (latitude, longitude, pres) float64 nan nan nan ... 35.05 35.05
    temp       (latitude, longitude, pres) float64 nan nan nan ... 3.934 3.933
    gamma_n    (latitude, longitude, pres) float64 nan nan nan ... 27.98 27.98
    SA         (latitude, longitude, pres) float64 nan nan nan ... 35.22 35.22
    ...         ...
    Kt         (latitude, longitude, pres) float64 nan nan ... 2.579e-07
    KtTz       (latitude, longitude, pres) float64 nan nan ... 1.215e-09
    Tz~        (latitude, longitude, pres) float64 nan nan nan ... nan nan nan
    chi~       (latitude, longitude, pres) float64 nan nan nan ... nan nan nan
    Kt~        (latitude, longitude, pres) float64 nan nan nan ... nan nan nan
    KtTz~      (latitude, longitude, pres) float64 nan nan nan ... nan nan nan
Attributes: (12/13)
    Conventions:           CF-1.6
    netcdf_version:        4
    project:               North Atlantic Tracer Release Experiment (NATRE)
    expocode:              32OC250_4
    cast_number:           3.0
    title:                 Microstructure profiler data from the ship Oceanus...
    ...                    ...
    latitude:              27.533166666666666
    longitude:             -30.723333333333333
    chief_scientist:       Raymond W. Schmitt
    data_originator:       Polzin
    institution:           WHOI
    data_assembly_center:  CCHDO
from xarray.tests import raise_if_dask_computes

bins = ed.sections.choose_bins(natre.gamma_n, depth_range=np.arange(150, 2001, 150))
with raise_if_dask_computes():
    micro150 = ed.sections.bin_average_vertical(
        natre.stack(
            {"cast": ["latitude", "longitude"]}, create_index=False
        ).assign_coords(cast=("cast", np.arange(100), {"cf_role": "profile_id"})),
        "neutral_density",
        bins,
        blocksize=20,
    )
2022-06-01 15:27:01,089 - distributed.nanny - WARNING - Restarting worker
2022-06-01 15:27:01,354 - distributed.nanny - WARNING - Restarting worker
micro150.load(client=client)
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/statsmodels/compat/pandas.py:65: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
  from pandas import Int64Index as NumericIndex
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/statsmodels/compat/pandas.py:65: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
  from pandas import Int64Index as NumericIndex
<xarray.Dataset>
Dimensions:           (bound: 2, gamma_n: 12, cast: 100)
Coordinates:
  * bound             (bound) <U5 'lower' 'upper'
  * gamma_n           (gamma_n) float64 26.61 26.91 27.12 ... 27.9 27.94 27.96
    num_obs           (gamma_n) int64 27405 28834 27525 ... 30729 29782 21047
    pres              (gamma_n) float64 225.0 375.1 ... 1.746e+03 1.878e+03
    gamma_n_bounds    (bound, gamma_n) float64 26.43 26.78 27.03 ... 27.95 27.97
    latitude          (cast) float64 27.5 27.5 27.5 27.5 ... 23.9 23.9 23.9 23.9
    longitude         (cast) float64 -30.7 -30.3 -29.8 ... -27.6 -27.2 -26.8
    time              (cast) datetime64[ns] 1992-03-28T15:28:59.999996928 ......
  * cast              (cast) int64 0 1 2 3 4 5 6 7 8 ... 92 93 94 95 96 97 98 99
Data variables: (12/59)
    chi               (gamma_n) float64 1.538e-08 4.209e-09 ... 2.571e-10
    eps               (gamma_n) float64 9.01e-10 5.513e-10 ... 1.245e-10
    KtTz              (gamma_n) float64 1.978e-07 1.337e-07 ... 2.619e-08
    KtTz~             (gamma_n) float64 2.604e-07 1.44e-07 ... 5.607e-08
    chib2             (gamma_n) float64 7.689e-09 2.105e-09 ... 1.285e-10
    hm                (gamma_n) float64 146.5 152.9 146.5 ... 165.6 159.5 114.6
    ...                ...
    δKtTzTz           (gamma_n) float64 1.479e-09 3.49e-10 ... 1.771e-11
    δKtTz~Tz          (gamma_n) float64 6.41e-10 3.147e-10 ... 1.952e-11
    δwTTz             (gamma_n) float64 nan nan nan nan nan ... nan nan nan nan
    δresidual         (gamma_n) float64 7.685e-10 2.61e-11 ... 9.237e-12
    δresidual_chi     (gamma_n) float64 2.518e-10 -4.472e-11 ... -2.887e-12
    δpres             (gamma_n) float64 73.26 76.45 73.26 ... 82.8 79.77 57.32
Attributes: (12/14)
    Conventions:           CF-1.6
    netcdf_version:        4
    project:               North Atlantic Tracer Release Experiment (NATRE)
    expocode:              32OC250_4
    cast_number:           3.0
    title:                 Microstructure profiler data from the ship Oceanus...
    ...                    ...
    longitude:             -30.723333333333333
    chief_scientist:       Raymond W. Schmitt
    data_originator:       Polzin
    institution:           WHOI
    data_assembly_center:  CCHDO
    commit:                eddydiff: b55f9526739dd63a768df2f816bc8327bcc0d247...

Latest one with “sign stabilized” \(⟨K_T T_z⟩\) (purple). For some reason this one is sometimes larger than \(⟨χ⟩\) in the upper bins

ed.plot.debug_section_estimate(micro150, KρTz2=True)
_images/natre-examine_57_0.png
micro150.to_netcdf("../datasets/natre-density-binned.nc")

Older#

Notation:

  1. "KtTz" is “naive”

  2. "KtTz~" is estimated over 5m bins. I try to make the sign of the gradient somewhat “stable”.

  3. "wTTz" is estimated in density space, without accounting for sign of \(\widetilde{T_z}\)

Latest one with “sign stabilized” \(⟨K_T T_z⟩\). For some reason this one is larger than \(⟨χ⟩\) in the upper bins

ed.plot.debug_section_estimate(micro150, KρTz2=True)
_images/natre-examine_61_0.png

Bigger error bars with 100m bins. Averaging to 150m bins reduces the error.

dcpy.plots.fill_between_bounds(micro100, "KtTz", y="pres")
dcpy.plots.fill_between_bounds(micro150, "KtTz", y="pres")
plt.xscale("log")
_images/natre-examine_63_0.png
dcpy.plots.fill_between_bounds(micro150, "chib2", y="pres")
dcpy.plots.fill_between_bounds(micro150, "KρTz2", y="pres")
dcpy.plots.fill_between_bounds(micro150, "KtTzTz", y="pres")
# dcpy.plots.fill_between_bounds(micro, "wTTz", y="pres")
plt.xlim([1e-12, 1e-8])
plt.xscale("log")
_images/natre-examine_64_0.png

THe decrease in “naive” 0.5-m \(⟨K_T T_z⟩\) is because I was using \(T_z\) > 1e-3. This ends up throwing out a lot of points near the bottom.

dcpy.plots.fill_between_bounds(micro, "chib2", y="pres")
dcpy.plots.fill_between_bounds(micro, "KρTz2", y="pres")
dcpy.plots.fill_between_bounds(micro, "KtTzTz", y="pres")
# dcpy.plots.fill_between_bounds(micro, "wTTz", y="pres")
plt.xlim([1e-12, 1e-8])
plt.xscale("log")
_images/natre-examine_66_0.png

The black is estimated in density space without accounting for sign of \(\widetilde{T_z}\), so we recover χ

dcpy.plots.fill_between_bounds(micro, "chib2", y="pres")
dcpy.plots.fill_between_bounds(micro, "KρTz2", y="pres")
dcpy.plots.fill_between_bounds(micro, "KtTzTz", y="pres")
dcpy.plots.fill_between_bounds(micro, "wTTz", y="pres")
plt.xlim([1e-12, 1e-8])
plt.xscale("log")
_images/natre-examine_68_0.png

Different bootstrapping method#

Explicitly calculate degrees of freedom like in Ferrari & Polzin (2005)

# micro.pres.attrs["bounds"] = "pres_err"
for label, group in grouped:
    print(label)
    break
(27.74, 27.79]
ds = group.unstack()
da = ds.KtTz

dp = 0.5
corrscales = [[slice(0, 800), 5], [slice(800, None), 10]]
dof = ed.sections.get_dof(da, dp=0.5, corrscales=corrscales)
dof
1011
from arch.bootstrap import IIDBootstrap

array = da.data.reshape(-1)
print("old...", ed.sections.compute_bootstrapped_mean_ci(array, 20, clean=False))
array = array[~np.isnan(array)]

plt.plot(array, marker="x")
absarray = array  # np.abs(array)
thresh = np.mean(absarray) + 30 * np.std(absarray)
array = np.where(np.abs(array) < thresh, array, np.nan)
array = array[~np.isnan(array)]

plt.plot(array, marker="x")

mean = np.mean(array)
stderr = np.mean(IIDBootstrap(array[~np.isnan(array)]).apply(np.std)) / np.sqrt(dof)
print("new...", [mean - stderr, mean, mean + stderr])
ed.sections.compute_bootstrapped_mean_ci(array, 20, clean=False)
old... [4.50550872e-08 6.73010347e-08 1.42173197e-07]
new... [1.784538731438142e-08, 4.454195776804866e-08, 7.12385282217159e-08]
array([3.21298005e-08, 4.45419578e-08, 6.02803924e-08])
_images/natre-examine_72_2.png

Partitioning \(θ_z\)#

def bin_coarsen(ds, func="mean"):
    return getattr(
        ds.coarsen(latitude=10, longitude=10, pres=300, boundary="trim"), func
    )()
natre = ed.natre.read_natre(load=True).sel(pres=slice(2000))
natre_binned = xr.load_dataset("../datasets/natre-density-binned.nc")
natre_coarse = natre.where(natre.chi_masked.notnull()).pipe(bin_coarsen).squeeze()
f, ax = plt.subplots(1, 4, sharey=False)

for Tzmin in [1e-3, 2e-3, 5e-3]:
    mask = (np.abs(natre.Tz) > Tzmin) & (natre.chi < 1e-6)
    chib2 = (natre.chi / 2).where(mask)
    Kt = chib2 / natre.Tz.where(mask) ** 2
    KtTz = chib2 / natre.Tz.where(mask)

    bin_coarsen(chib2).squeeze().cf.plot(
        y="pres", label=str(Tzmin), ax=ax[0], xscale="log"
    )

    bin_coarsen(Kt).squeeze().cf.plot(
        y="pres", label=str(Tzmin), ax=ax[1], xscale="log"
    )

    bin_coarsen(KtTz).squeeze().cf.plot(
        y="pres", label=str(Tzmin), ax=ax[2], xscale="log"
    )

    (bin_coarsen(KtTz, func="std") / np.sqrt(100 * 10)).squeeze().cf.plot(
        y="pres", label=str(Tzmin), ax=ax[3], xscale="log"
    )

    # np.log(1025 * 4000 * np.abs(KtTz)).plot.hist(
    #    bins=np.arange(-10, 10, 0.05), histtype="step", ax=ax[3], density=True
    # )

(natre_binned.chib2).cf.plot(y="pres", ax=ax[0], color="k")
(natre_binned.KtTzTz / natre_binned.dTdz_m).cf.plot(y="pres", ax=ax[2])
(natre_binned.chib2 / natre_binned.dTdz_m).cf.plot(y="pres", ax=ax[2], color="k")
(natre_binned.Krho_m * natre_binned.dTdz_m).cf.plot(y="pres", ax=ax[2], color="c")

(natre_binned.δKtTzTz / natre_binned.dTdz_m).cf.plot(y="pres", ax=ax[3], color="k")
(natre_binned.δKρTz2 / natre_binned.dTdz_m).cf.plot(y="pres", ax=ax[3], color="k")
ax[1].legend()
dcpy.plots.clean_axes(ax)
f.set_size_inches((9, 4))
_images/natre-examine_76_0.png
Kt = natre.chi / 2 / natre.Tz.where(np.abs(natre.Tz) > 5e-3) ** 2
natre_coarse.Kt.plot(y="pres")
(natre_binned.Krho_m).plot(y="pres", color="k")
(natre_binned.Kt_m / natre_binned.dTdz_m).cf.plot(y="pres", xscale="log")
bin_coarsen(Kt).plot(y="pres")
# Kt.coarsen({"latitude": 10, "longitude": 10, "pres_3.1": 30}, boundary="trim").mean().plot(y="pres_3.1")
[<matplotlib.lines.Line2D at 0x16920d910>]
_images/natre-examine_77_1.png
(bin_coarsen(Kt) * natre_coarse.temp.differentiate("pres") ** 2).squeeze().cf.plot()
(natre_coarse.Kt * natre_coarse.temp.differentiate("pres") ** 2).cf.plot(xscale="log")
(natre_coarse.chi / 2).cf.plot(xscale="log")
natre_binned.chib2.plot(y="pres")
[<matplotlib.lines.Line2D at 0x168e38100>]
_images/natre-examine_78_1.png
dTdz.coarsen(
    {"pres_3.1": 61, "latitude": 10, "longitude": 10}, boundary="trim"
).mean().plot()
natre_binned.dTdz_m.plot(x="pres", yscale="log")
[<matplotlib.lines.Line2D at 0x178c62370>]
_images/natre-examine_79_1.png
natre_binned.dTdz_m.plot(y="pres")
(-1 * natre_coarse.temp.differentiate("pres")).plot(y="pres")
[<matplotlib.lines.Line2D at 0x16a09be20>]
_images/natre-examine_80_1.png

Compare neutral density bins to FP2005#

plt.figure()
fp = xr.DataArray(
    [26.59, 26.95, 27.22, 27.46, 27.64, 27.77, 27.85, 27.91, 27.95],
    dims=("depth",),
    coords={"depth": [200, 400, 600, 800, 1000, 1200, 1400, 1600, 1800]},
)
fp.plot(marker="x", label="FP2005")
micro.gamma_n.plot(x="pres", marker=".", label="mine")
# plt.plot(np.arange(150, 2001, 100), bins, marker='^')
plt.gca().set_xticks(fp.depth)
plt.legend()
<matplotlib.legend.Legend at 0x7f2cd6b69f70>
_images/natre-examine_82_1.png

Better \(\widetilde{K_T T_z}\) estimate#

It appears that preserving the sign of \(wT\) is important. After a lot of averaging I do get nearly all positive values but the error bars are big!

  • Yeah the sign for \(T_m + T_e\) can be negative and opposite to \(T_z^m\) Indeed that is important!

  • If I sort the \(Θ\) then I get rid of the stable inversions that are important.

  • ~No matter what I do, I basically recover χ :(.~ I recover χ if I don’t pay attention to this fact.

  • I think this explains why \(\widetilde{K_t T_z}\) with sign preserved is closer to \(Γ ⟨ε⟩/N_m^2 T_z^m\). The cancellations are important! Using absolute values kills the cancellations…

  • Might be better to stabilize the sign somehow; Logic here is that

    • the sign of \(\widetilde{T_z}\) is set by the intermediate scale (30-50m-ish);

    • but we want \(\widetilde{w_tθ_t} = \widetilde{χ}/\widetilde{T_z}/2\) over 5-m scales

    • So filter at 10m wavelength; choose sign;

    • Then determine \(\widetilde{T_z}\) magnitude over 5-m linear fits;

    • apply sign of first element in window (this is somewhat arbitrary)

natre = ed.natre.read_natre(load=True)
natre_binned = xr.load_dataset("../datasets/natre-density-binned.nc")
Tzmi = natre_binned.dTdz_m.reset_coords(drop=True).interp(gamma_n=natre.gamma_n)
Tmi = natre_binned.theta.reset_coords(drop=True).interp(gamma_n=natre.gamma_n)

Depth-space#

Here I estimate \(\widetilde{w^tθ^t} = \widetilde{χ} / 2 / \widetilde{T}_z\) over 5m scales.

  1. Using linear fit to \(T\) (or sorted \(T\)) ofr gradient;

  2. simple 5m average of χ

The core problem is that a simple averaging of \(\widetilde{T}_z\) basically yields \(T_z^m\); but it should be bigger for this strategy to work.

(natre.Tz - Tzmi).plot.line(row="latitude", col="longitude")
<xarray.plot.facetgrid.FacetGrid at 0x17cf48d60>
_images/natre-examine_89_1.png
import hvplot.xarray
import xfilter

loc = dict(latitude=3, longitude=8)

T = natre.CT.isel(loc).interpolate_na("pres").dropna("pres")
# Tm = xfilter.lowpass(T, coord="pres", freq=1 / 200, order=1)
Te = xfilter.lowpass(T, coord="pres", freq=1 / 20, order=1)

# (T.plot())
# (Te.plot())
# (Tmi.isel(loc).plot())
# xfilter.lowpass(
#    (T - Tmi.isel(loc)).interpolate_na("pres"), coord="pres", freq=1 / 20
# ).plot()
# (Tm.plot())
# ((T - Te).plot(ax=plt.gca().twinx()))

# (T.hvplot.line(x="pres") * Te.hvplot.line(x="pres"))
(T.differentiate("pres").hvplot.line()) * Te.differentiate("pres").hvplot.line()
clean.dTdz.isel(latitude=2, longitude=6).plot()
clean.dTdz.isel(latitude=9, longitude=8).plot()
natre_binned.dTdz_m.plot(x="pres", color="k")
[<matplotlib.lines.Line2D at 0x17c902df0>]
_images/natre-examine_91_1.png
(clean.dTdz).coarsen(
    {"latitude": 10, "longitude": 10, "pres_5.1": 20}, boundary="trim"
).mean().plot(yscale="log")
(natre_binned.dTdz_m).plot(x="pres")
[<matplotlib.lines.Line2D at 0x17757ff40>]
_images/natre-examine_92_1.png
clean = ed.sections.estimate_microscale_stirring_depth_space(
    natre, filter_len=10, segment_len=3
)

plt.figure()
clean["Tz~"].coarsen(pres=200, latitude=10, longitude=10, boundary="trim").mean().plot()
natre_binned.dTdz_m.plot(x="pres", color="k", yscale="log")

plt.figure()
natre.KtTz.where(np.abs(natre.Tz) > 1e-5).mean(["latitude", "longitude"]).plot(
    yscale="log", lw=0.5, label="simple"
)
(
    (
        natre.KtTz.where(np.abs(natre.Tz) > 1e-5)
        .coarsen(pres=200, latitude=10, longitude=10, boundary="trim")
        .mean()
    ).plot(yscale="log", label="simple coarsened")
)

(
    clean["KtTz~"]
    # .where(np.abs(clean["KtTz~"]) < 5e-6)
    .cf.coarsen({"Z": 200, "latitude": 10, "longitude": 10}, boundary="trim").mean()
    # .mean(["latitude", "longitude"])
    .plot(label="cleaner")
)

(natre_binned.KρTz2 / natre_binned.dTdz_m).plot(label="$K_ρT_z^m$", color="k", x="pres")
(natre_binned.chib2 / natre_binned.dTdz_m).plot(
    label="$χ/2/T_z^m$",
    color="k",
    x="pres",
    ylim=(1e-8, 1e-6),
)
# (micro.wTTz / micro.dTdz_m).plot(label="$K_tT_z$", color="cyan", x="pres")
plt.legend()
7 20
50243  values < 0
<matplotlib.legend.Legend at 0x177280df0>
_images/natre-examine_93_2.png _images/natre-examine_93_3.png
kwargs = dict(pres=100, latitude=10, longitude=10, boundary="trim")

clean.KtTz.mean(["latitude", "longitude"]).coarsen(
    {"pres_5.1": 20}, boundary="trim"
).mean().plot()

(
    ((natre.KtTz).where(np.abs(natre.Tz) > 1e-3).coarsen(**kwargs).mean()).plot(
        yscale="log"
    )
)

((natre.KrhoTz).coarsen(**kwargs).mean().plot(yscale="log", color="k"))
# np.abs(natre.chi / 2).mean(["latitude", "longitude"]).coarsen(**kwargs).mean().plot(
#    yscale="log"
# )
---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Input In [316], in <cell line: 3>()
      1 kwargs = dict(pres=100, latitude=10, longitude=10, boundary="trim")
----> 3 clean.KtTz.mean(["latitude", "longitude"]).coarsen(
      4     {"pres_5.1": 20}, boundary="trim"
      5 ).mean().plot()
      7 (
      8     (((natre.KtTz).where(np.abs(natre.Tz) > 1e-3).coarsen(**kwargs).mean())).plot(
      9         yscale="log"
     10     )
     11 )
     13 ((natre.KrhoTz).coarsen(**kwargs).mean().plot(yscale="log", color="k"))

NameError: name 'clean' is not defined
natre.eps.coarsen(**kwargs).mean().plot(yscale="log")
natre.eps_chi.where(np.abs(natre.Tz) > 1e-3).coarsen(**kwargs).mean().plot(yscale="log")
[<matplotlib.lines.Line2D at 0x16469f1f0>]
_images/natre-examine_95_1.png

Interestingly excluding small gradients excludes small χ, so the \(\widetilde{K_TT_z}\) estimate increases with more filtering!

kwargs = dict(bins=100, density=True, histtype="step")
np.log10(chi / 2 / dTdz).plot.hist(**kwargs)
np.log10(chi / 2 / dTdz.where(dTdz > 5e-3)).plot.hist(**kwargs);
_images/natre-examine_97_0.png

T-space#

This is not as easy.

  1. Choosing T bins is a real pain. Because we are covering a giant region; the spread in cast-mean T is large. So picking bins is a little challening.

  2. By averaging over all casts I end up with \(⟨χ⟩\), and without accounting for sign of \(Δz\) (separation between isotherms), I end up with \(T_z^m\) in the denomoinator.

natre = ed.natre.read_natre(load=True)
clean = ed.sections.estimate_microscale_stirring_depth_space(
    natre.stack({"cast": ("latitude", "longitude")}).drop("cast"),  # .isel(cast=36),
    filter_len=10,
    segment_len=5,
)
clean
11 20
/Users/dcherian/mambaforge/envs/eddydiff/lib/python3.9/site-packages/numba/np/ufunc/gufunc.py:151: RuntimeWarning: invalid value encountered in get_gradient_sign
  return self.ufunc(*args, **kwargs)
0  values == 0
33882
<xarray.Dataset>
Dimensions:  (pres: 3981, cast: 100)
Coordinates:
  * pres     (pres) float64 10.0 10.5 11.0 11.5 ... 1.999e+03 2e+03 2e+03
  * cast     (cast) int64 0 1 2 3 4 5 6 7 8 9 ... 90 91 92 93 94 95 96 97 98 99
    window   (pres) float64 -0.0 -0.5 -1.0 -1.5 -2.0 ... nan nan nan nan nan
Data variables:
    Tz~      (pres, cast) float64 nan nan nan nan nan ... nan nan nan nan nan
    chi~     (pres, cast) float64 nan nan nan nan nan ... nan nan nan nan nan
    Kt~      (pres, cast) float64 nan nan nan nan nan ... nan nan nan nan nan
    KtTz~    (pres, cast) float64 nan nan nan nan nan ... nan nan nan nan nan
    T~       (cast, pres) float64 nan nan nan nan nan ... nan nan nan nan nan
natre = natre.update(clean)
from xarray.tests import raise_if_dask_computes

bins = ed.sections.choose_bins(natre.gamma_n, depth_range=np.arange(150, 2001, 100))
with raise_if_dask_computes():
    micro = ed.sections.bin_average_vertical(
        natre.cf.stack(
            {"cast": ("latitude", "longitude")}, create_index=False
        ).assign_coords({"cast": ("cast", np.arange(100), {"cf_role": "profile_id"})}),
        "neutral_density",
        bins,
        blocksize=10,
    )
# micro.pres.attrs["bounds"] = "pres_err"
micro
<xarray.Dataset>
Dimensions:         (bound: 2, gamma_n: 18, cast: 100)
Coordinates:
  * bound           (bound) <U5 'lower' 'upper'
  * gamma_n         (gamma_n) float64 26.56 26.78 26.95 ... 27.93 27.95 27.96
    num_obs         (gamma_n) int64 19358 19463 20127 ... 20191 11530 22215
    pres            (gamma_n) float64 201.4 300.6 401.1 ... 1.795e+03 1.877e+03
    gamma_n_bounds  (bound, gamma_n) float64 26.43 26.69 26.87 ... 27.95 27.97
    latitude        (cast) float64 27.5 27.5 27.5 27.5 ... 23.9 23.9 23.9 23.9
    longitude       (cast) float64 -30.7 -30.3 -29.8 -29.4 ... -27.6 -27.2 -26.8
    time            (cast) datetime64[ns] 1992-03-28T15:28:59.999996928 ... 1...
  * cast            (cast) int64 0 1 2 3 4 5 6 7 8 ... 92 93 94 95 96 97 98 99
Data variables: (12/54)
    chi             (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    eps             (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    KtTz            (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    KtTz~           (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    hm              (gamma_n) float64 98.77 98.73 101.8 ... 101.3 58.27 114.8
    wT              (gamma_n) float64 nan nan nan nan nan ... nan nan nan nan
    ...              ...
    δKtTz~Tz        (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    δwTTz           (gamma_n) float64 nan nan nan nan nan ... nan nan nan nan
    δresidual       (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    δpres           (gamma_n) float64 49.38 49.37 50.92 ... 50.63 29.14 57.4
    chib2           (gamma_n) float64 dask.array<chunksize=(1,), meta=np.ndarray>
    chib2_err       (gamma_n, bound) float64 dask.array<chunksize=(1, 2), meta=np.ndarray>
Attributes: (12/14)
    Conventions:           CF-1.6
    netcdf_version:        4
    project:               North Atlantic Tracer Release Experiment (NATRE)
    expocode:              32OC250_4
    cast_number:           3.0
    title:                 Microstructure profiler data from the ship Oceanus...
    ...                    ...
    longitude:             -30.723333333333333
    chief_scientist:       Raymond W. Schmitt
    data_originator:       Polzin
    institution:           WHOI
    data_assembly_center:  CCHDO
    commit:                eddydiff: 3c5e37fac5bc4c53f9567addc4c058ef415a63b4...
micro.load()
<xarray.Dataset>
Dimensions:         (bound: 2, gamma_n: 18, cast: 100)
Coordinates:
  * bound           (bound) <U5 'lower' 'upper'
  * gamma_n         (gamma_n) float64 26.56 26.78 26.95 ... 27.93 27.95 27.96
    num_obs         (gamma_n) int64 19358 19463 20127 ... 20191 11530 22215
    pres            (gamma_n) float64 201.4 300.6 401.1 ... 1.795e+03 1.877e+03
    gamma_n_bounds  (bound, gamma_n) float64 26.43 26.69 26.87 ... 27.95 27.97
    latitude        (cast) float64 27.5 27.5 27.5 27.5 ... 23.9 23.9 23.9 23.9
    longitude       (cast) float64 -30.7 -30.3 -29.8 -29.4 ... -27.6 -27.2 -26.8
    time            (cast) datetime64[ns] 1992-03-28T15:28:59.999996928 ... 1...
  * cast            (cast) int64 0 1 2 3 4 5 6 7 8 ... 92 93 94 95 96 97 98 99
Data variables: (12/54)
    chi             (gamma_n) float64 1.664e-08 5.465e-09 ... 2.303e-10
    eps             (gamma_n) float64 9.838e-10 5.747e-10 ... 1.235e-10
    KtTz            (gamma_n) float64 2.515e-07 1.472e-07 ... 3.07e-08 2.191e-08
    KtTz~           (gamma_n) float64 2.887e-07 1.79e-07 ... 5.764e-08 5.291e-08
    hm              (gamma_n) float64 98.77 98.73 101.8 ... 101.3 58.27 114.8
    wT              (gamma_n) float64 nan nan nan nan nan ... nan nan nan nan
    ...              ...
    δKtTz~Tz        (gamma_n) float64 7.741e-10 4.478e-10 ... 2.664e-11 1.7e-11
    δwTTz           (gamma_n) float64 nan nan nan nan nan ... nan nan nan nan
    δresidual       (gamma_n) float64 7.387e-10 5.777e-11 ... 6.443e-12
    δpres           (gamma_n) float64 49.38 49.37 50.92 ... 50.63 29.14 57.4
    chib2           (gamma_n) float64 8.32e-09 2.733e-09 ... 1.428e-10 1.151e-10
    chib2_err       (gamma_n, bound) float64 7.713e-09 8.926e-09 ... 1.229e-10
Attributes: (12/14)
    Conventions:           CF-1.6
    netcdf_version:        4
    project:               North Atlantic Tracer Release Experiment (NATRE)
    expocode:              32OC250_4
    cast_number:           3.0
    title:                 Microstructure profiler data from the ship Oceanus...
    ...                    ...
    longitude:             -30.723333333333333
    chief_scientist:       Raymond W. Schmitt
    data_originator:       Polzin
    institution:           WHOI
    data_assembly_center:  CCHDO
    commit:                eddydiff: 3c5e37fac5bc4c53f9567addc4c058ef415a63b4...

Using fancy sign detection#

def notnull(data):
    flat = data.data.ravel()
    return flat[~np.isnan(flat)]
mask = np.abs(clean["KtTz~"]) > 1e-5
clean["chi~"].where(mask).pipe(notnull)
array([2.29729636e-05, 4.99936064e-06, 8.38180596e-08, 4.42592055e-08,
       2.45233000e-08, 6.14768901e-08, 5.21607218e-07, 2.45605935e-07,
       9.83617336e-08, 2.11507639e-07, 1.67454668e-08])
clean["Tz~"].where(np.abs(clean["KtTz~"]) > 1e-5).pipe(notnull)
array([ 0.14074649,  0.13337417,  0.00319876,  0.00175999,  0.00105803,
        0.00300432,  0.00974087,  0.00935505, -0.00478368,  0.00899659,
        0.00078514])
clean["Tz~"].isel(cast=16).plot(marker="x")
[<matplotlib.lines.Line2D at 0x191bccb80>]
_images/natre-examine_108_1.png
clean["KtTz~"].isel(cast=26).ffill("pres").plot(yscale="linear")
[<matplotlib.lines.Line2D at 0x193737850>]
_images/natre-examine_109_1.png
ed.plot.debug_section_estimate(micro)
_images/natre-examine_110_0.png
dcpy.plots.fill_between_bounds(micro, "chib2", y="pres")
dcpy.plots.fill_between_bounds(micro, "KρTz2", y="pres")
dcpy.plots.fill_between_bounds(micro, "KtTz~Tz", y="pres")
# dcpy.plots.fill_between_bounds(micro, "KtTzTz", y="pres")
# dcpy.plots.fill_between_bounds(micro, "wTTz", y="pres")
plt.xlim([1e-11, 1e-8])
plt.xscale("log")
_images/natre-examine_111_0.png

Single bin#

from xarray.tests import raise_if_dask_computes

bins = ed.sections.choose_bins(natre.gamma_n, depth_range=np.arange(150, 2001, 100))
with raise_if_dask_computes():
    grouped = ed.sections.bin_average_vertical(
        natre.cf.stack(
            {"cast": ("latitude", "longitude")}, create_index=False
        ).assign_coords({"cast": ("cast", np.arange(100), {"cf_role": "profile_id"})}),
        "neutral_density",
        bins[10:],
        blocksize=10,
        return_group=True,
    )
# micro.pres.attrs["bounds"] = "pres_err"
for label, group in grouped:
    break
Tbins, sortT = ed.sections.estimate_microscale_stirring(
    group.unstack(), dz=7, debug=True
)
dT = np.diff(Tbins).mean()
sort = dcpy.oceans.thorpesort(sortT, by="gamma_n", core_dim="pres")
sort.CT.plot.line(hue="cast", add_legend=False);
_images/natre-examine_114_0.png
sort.CT.isel(cast=[20, 30, 60]).cf.dropna("Z").plot(hue="cast")
dcpy.plots.liney(Tbins)
_images/natre-examine_115_0.png
CT = profile.CT.data
profiles = natre.copy(deep=True)
profiles["CT"] = profiles.CT.interpolate_na("pres")
profiles["Tzfilt"] = -1 * xfilter.lowpass(
    profiles.CT, coord="pres", freq=1 / 10
).differentiate("pres").bfill("pres")
stacked = profiles.stack({"cast": ("latitude", "longitude")}, create_index=False)
profile = stacked.interpolate_na("pres")  # .sel(pres=slice(1120, 1250))


# get_gradient_sign(
#    profile.pres.data, profile.CT.data, profile.Tzfilt.data, 5, 10, a.data
# )
# profile.CT.plot()
Tzsign.plot(ax=plt.gca().twinx())
<matplotlib.collections.QuadMesh at 0x185d14790>
_images/natre-examine_119_1.png
def clean_idx(idx, CT, ΔP, window):

    P = CT["pres"].data
    CT = CT.data

    def _clean(idx, P, ΔP, window):
        newidx = [idx[0]]
        for left, right in zip(idx[:-1], idx[1:]):
            Pl, Pr = P[left], P[right]
            # print(np.abs(Pr - Pl))
            if np.abs(Pr - Pl) > ΔP:
                if ((right - left) % window) < window // 2:
                    right += 1
                newidx.append(right)

        if newidx[-1] != idx[-1]:
            # Last index was dropped
            # instead to preserve all data;
            # we drop the previous one and re-add the last one
            newidx.pop()
            newidx.append(idx[-1])
        print(idx, newidx)
        return newidx

    idx = _clean(idx, P, ΔP, window=window)
    idx = _clean(idx, CT, 1e-4 * ΔP, window=1e10)
    return idx


def smooth_gradient(profile, chi, ΔP, window, debug=False):

    dp = 0.5
    smooth = xfilter.lowpass(
        profile, coord="pres", freq=1 / int(ΔP // dp), num_discard=5
    )
    Tzsmooth = smooth.differentiate("pres")
    # array = profile.data
    if debug:
        f, ax = plt.subplots(4, 1, sharex=True, constrained_layout=True)
        profile.plot(ax=ax[0])
        smooth.plot(ax=ax[0])

    (idx,) = np.nonzero(np.abs(np.sign(Tzsmooth).diff("pres").data) > 0)
    idx = np.concatenate([[0], idx, [profile.sizes["pres"] - 1]])
    if debug:
        dcpy.plots.linex(Tzsmooth["pres"][idx], color="k", lw=0.5, ax=ax[0])

    idx = clean_idx(idx, profile, ΔP=ΔP, window=window)
    if debug:
        dcpy.plots.linex(Tzsmooth["pres"][idx], color="b", ax=ax[0])

    Tzs = []
    for subprof in np.split(profile, idx[1:-1]):
        # -1 for temperature slope
        sign = np.sign(subprof[0] - subprof[-1]).data
        coarse = subprof.coarsen(pres=window, boundary="pad").construct(
            pres=("p", "window")
        )
        sorty = dcpy.oceans.thorpesort(
            coarse, coarse, core_dim="window", ascending=False
        )
        sorty = sorty.where(sorty.count("window") > 2)
        Tzsubprof = sign * np.abs(
            sorty.polyfit("window", deg=1)
            .sel(degree=1, drop=True)
            .polyfit_coefficients.data
        )
        Tzs.append(np.repeat(Tzsubprof, window)[..., : subprof.shape[-1]])

    dTdz = profile.copy(data=np.concatenate(Tzs)).ffill("pres")
    if debug:
        dTdz.plot(ax=ax[1])
        (-1 * Tzsmooth).plot(ax=ax[1])

        (chi).plot(yscale="log", ylim=(1e-12, None), ax=ax[2])
        ax2 = ax[2].twinx()
        (dTdz).plot(ax=ax2, color="k")
        (dTdz.where(np.abs(dTdz) < 2e-4)).plot(color="b", marker="o", ls="none", ax=ax2)
        dcpy.plots.set_axes_color(ax2, "k")

        (1025 * 4000 * -chi / dTdz).plot(ax=ax[3], yscale="symlog")
        (1025 * 4000 * -chi / -Tzsmooth).plot(ax=ax[3])
        (1025 * 4000 * -chi / dTdz.where(np.abs(dTdz) < 2e-4)).plot(
            color="b", marker="o", ls="none", ax=ax[3]
        )

        dcpy.plots.clean_axes(np.atleast_2d(ax).T)
        f.set_size_inches((10, 7))

    return dTdz
    # dcpy.plots.liney(Tbins)


# NATRE 36, 56, 60 are good to test with
profile = sort.isel(cast=36).cf.dropna("Z")

dTdz = smooth_gradient(profile.CT, profile.chi, ΔP=5, window=20, debug=True);
[  0   6  11  37  59  76 104 119 124 135 148 153 187] [0, 38, 60, 76, 105, 119, 135, 148, 187]
[0, 38, 60, 76, 105, 119, 135, 148, 187] [0, 39, 61, 77, 106, 120, 136, 149, 187]
_images/natre-examine_120_1.png
Tz_m = sort.CT.polyfit("pres", deg=1).sel(degree=1, drop=True).polyfit_coefficients
Tz_m.plot()
[<matplotlib.lines.Line2D at 0x180d0ceb0>]
_images/natre-examine_121_1.png
Tz = -1 * xfilter.lowpass(sort.CT, coord="pres", freq=1 / 10).differentiate("pres")
KtTz = sort.chi / Tz
KtTz = KtTz.where(np.abs(KtTz) < 300 / 1025 / 4000)
(KtTz.mean("pres") * Tz_m).plot()
(sort.chi.mean("pres") / 2).plot()
[<matplotlib.lines.Line2D at 0x1814e42b0>]
_images/natre-examine_123_1.png

Try ε#

natre.eps_chi.where(np.abs(natre.Tz) > 1e-10).mean(["latitude", "longitude"]).plot(
    yscale="log"
)
natre.eps_chi.where(np.abs(natre.Tz) > 1e-2).mean(["latitude", "longitude"]).plot(
    yscale="log"
)
natre.eps.mean(["latitude", "longitude"]).plot(yscale="log")
[<matplotlib.lines.Line2D at 0x170f51220>]
_images/natre-examine_125_1.png

Sensitivity of χpod estimate to mean K estimate#

Various estimates of the variance produced by turbulent stirring of mean vertical gradient

  1. Ferrai & Polzin (2005): \(Γ ⟨ε⟩/N_m^2 ∂_zθ_m^2\)

  2. χpod estimate: \(⟨K_T θ_z⟩ ∂_zθ_m\)

  3. (1) with \(ε_χ\) : \(Γ ⟨ε_χ⟩/N_m^2 ∂_zθ_m^2\)

  4. (2) with \(K_ρ\) instead of \(K_T\): \(⟨K_ρ θ_z⟩ ∂_zθ_m^2\)

(chidens.Krho_m * chidens.dTdz_m**2).cf.plot.step(
    color="k", lw=2, label="$Γ ⟨ε⟩/N_m^2  ∂_zθ_m^2$"
)

(chidens.KtTz * chidens.dTdz_m).cf.plot.step(
    color="k", marker="x", ls="none", label="$⟨K_T θ_z⟩ ∂_zθ_m$"
)

(0.2 * chidens.eps_chi / chidens.N2_m * chidens.dTdz_m**2).cf.plot.step(
    color="C1", lw=2, label="$Γ ⟨ε_χ⟩/N_m^2  ∂_zθ_m^2$"
)

(chidens.KrhoTz * chidens.dTdz_m).cf.plot.step(
    lw=2, label="$⟨K_ρ θ_z⟩  ∂_zθ_m$", color="C2"
)

plt.grid(True, which="both", lw=0.5)

plt.legend()
plt.xscale("log")
plt.xlabel("Variance production or dissipation [°C²/s]")
plt.gcf().set_size_inches((4, 5))
_images/natre-examine_128_0.png
f, ax = plt.subplots(1, 2, sharey=True, sharex=True)
micro.Krho.cf.plot.step(y="pres", xscale="log", ax=ax[0])
micro.Kt.cf.plot.step(y="pres", ax=ax[0])
micro.Krho_m.cf.plot.step(y="pres", ax=ax[1])
micro.Kt_m.cf.plot.step(y="pres", ax=ax[1])
ax[0].legend(["$⟨K_ρ⟩$", "$⟨K_T⟩$"])
ax[1].legend(["$K_ρ^m = Γ ⟨ε⟩/N_m²$", "$K_T^m = ⟨χ⟩/2θ_m^2$"])
[aa.grid(True, which="both", lw=0.5) for aa in ax]
plt.gcf().set_size_inches((6, 5))
_images/natre-examine_129_0.png